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In this paper, we apply the probabilistic model checker PRISM to the analysis of a biological sys- 
tem - the Platelet-Derived Growth Factor (PDGF) signaling pathway, demonstrating in detail how 
this pathway can be analyzed in PRISM. We show that quantitative verification can yield a better 
understanding of the PDGF signaling pathway. 

1 Introduction 

Biological systems consist of separated components, which interact to influence each other and there- 
fore the whole system's behavior. The field of systems biology aims to understand such complex in- 
teractions from a systematic view. Due to the similarity between biological systems and complex dis- 
tributed/reactive systems studied in computer science 1301 , modeling and analyzing techniques developed 
in the field of formal methods can be applied to biological systems as well 0. Due to efficient verifica- 
tion techniques, formal methods can analyze large systems exhibiting complex behaviors - this process 
is typically supported by automatic computer tools. Clearly, this gives formal methods an advantage, 
as in silico experiments are much easier to perform than in vitro experiments for the aim of analyzing 
and understanding biological systems. During the last decade, there has been a rapid and successful de- 
velopment in applying formal methods to systems biology - new formalisms are developed for systems 
biology to create models for biological phenomena, new algorithms and tools are specially designed and 
tailored for the analysis of such models (e.g., see |[T0ll3Tl[TTl0 . 

In this paper, we explore the usage of model checkingfor biological systems. Model checking is 
referred to as the automatic process of checking whether a system model satisfies a given specification 
(expressed as a temporal logic formula), by exhaustively exploring all possible executions of the sys- 
tem. This differs from simulation-based techniques, which only study a subset of the executions. More 
specifically, we focus on the probabilistic model checking approach, first introduced by Hart, Shark 
and Pnueli ifTTl . as biological systems usually have complicated stochastic behaviors. This technique 
is well-established and widely used for ascertaining the correctness of real-life systems, including dis- 
tributed systems and communication protocols. In probabilistic model checking, systems are normally 
represented by Markov chains or Markov decision processes. Properties of the models are expressed 
in quantitative extensions of temporal logics. Probabilistic verification has gained notable success in 
analyzing probabilistic systems including biological signaling pathways (e.g., see Il25ll26l ). 
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We use the probabilistic model checker PRISM [26] to yield a better understanding of the PDGF 
(Platelet-Derived Growth Factor) signaling pathway. PDGF, described approximately 30 years ago as 
a major mitogenic component of whole blood ||37l , is a growth factor that regulates cell growth and 
division. It promotes angiogenesis and also preserves vascular integrity through the recruitment of per- 
icytes to endothelial tubes. Clinical studies reveal that aberrant expression of PDGF and its receptor 
is often associated with a variety of disorders such as atherosclerosis, fibroproliferative diseases and 
most importantly, neoplasia B71 . Deregulation of the PDGF signaling pathway plays a critical role in 
the development of many types of human diseases such as gastrointestinal stromal tumor and hypere- 
osinophilic syndrome lPT2l [Tl |3"71 151 l28l. Based on intensive literature review, we have built the PDGF 
signal transduction model in ODE (Ordinary Differential Equation) format. The essential part of the 
PDGF signaling pathway contains the coupling of PDGF ligand to its receptor PDGFR, the negative regu- 
latory mechanism on PDGFR and the activation of two main downstream signaling pathways, i.e., MAPK 
(Mitogen-Activated Protein Kinase) and PI3K/Akt pathways. In addition, there also exist positive and 
negative crosstalk interactions between different downstream signaling pathways (more details on the 
PDGF signaling pathway can be found in Section[3]>. In our study, there are three main goals: (1) ana- 
lyze the dynamics of PDGF induced signaling, (2) analyze the influence of the crosstalk reactions and 
(3) analyze the importance of individual reaction on downstream signaling molecules. The first two 
can be used to check whether the constructed signaling pathway is consistent with respect to biological 
data, while the last one can lead us to some prediction. We have achieved these goals by quantitative 
verification using PRISM. 

Related work. PRISM has been used for analyzing a variety of biological systems. In |fl8l , PRISM is used 
to analyze the FGF (Fibroblast Growth Factor) signaling pathway. Although only a model corresponding 
to a single instance of the pathway is built, it is still rich enough to explain the roles of the components 
in the pathway and how they interact. In ll24il . PRISM is used to study the MAPK cascade. The authors 
explain how the biological pathway can be modeled in PRISM and how this enables the analysis of a 
rich set of quantitative properties. 

Jha et al. [20] present the first algorithm for performing statistical model checking using Bayesian 
sequential hypothesis testing and test the performance of the algorithm on the FGF signaling pathway 
and several others. Schwarick and Heiner [32] give an interval decision diagram based approach to 
CSL model checking. They aim for efficient verification of biochemical network models and apply their 
approach to the MAPK cascade. Draabik et al. [9] apply the technique modular verification to analyzing 
the lac operon regulation. 

Apart from probabilistic (or stochastic) model checking, other formal techniques have been used for 
studying biological systems in recent years as well. For instance, the process algebra bio-PEPA is used to 
model and analyze the NF-kappaB pathway ; Petri nets are used for quantitative predictive modeling 
of C.elegans vulval development 10; and the rule-based modeling language Kappa is used to model 
epigenetic information maintenance [23]. More related work, especially case studies of applying formal 
methods in systems biology, can be found in [31 , 11]. 

Outline of the paper. In Section |2j we give an overview of probabilistic model checking and the tool 
PRISM. Section [3] describes the PDGF signaling pathway. In Section|4| we build a model in PRISM for 
the PDGF signaling pathway and describe several properties of the model that we are interested in. Our 
verification results are given in Section [5] Finally, we draw conclusion of this paper and discuss some 
future work in Section [6] 
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2 Probabilistic Model Checking and PRISM 

We briefly introduce probabilistic verification and the probabilistic model checker - PRISM [26]. 
2.1 CTMC and CSL 

Probabilistic model checking is a valiant of model checking, which aims at analyzing the correctness 
of finite state systems with focus on quantitative aspects. Model checking of a system requires two 
inputs: a formal description of the system, which is usually given in a high level modeling formalism 
(e.g., Petri nets or process algebra) and a specification of the system properties, which is usually given 
as temporal logic (e.g., CTL or LTL) formulas. After accepting the two inputs, a model checking tool 
then can verify whether the system satisfies the desired properties and give counter-examples if the 
system does not satisfy a certain property, by exploring all possible behaviors of the system exhaustively. 
As the word "probabilistic" indicates, probabilistic model checking focuses on systems with stochastic 
behaviors. Instead of asking the model checker "will the molecule become active in the end?", we can 
ask "what is the probability of the molecule being active at the steady state?" or "what is the probability 
of the molecule being active at time instant ??". In probabilistic model checking, systems are normally 
represented by Markov chains or Markov decision processes. In this paper, we use continuous-time 
Markov chains (CTMCs), to build the signaling pathway models. 

A CTMC can model both (continuous) real time and probabilistic choice by assigning rates at tran- 
sitions between states. The formal definition of a CTMC is given as follows. 

Definition 2.1 Let M>o denote the set of non-negative reals andAP be a fixed finite set of atomic propo- 
sitions. A CTMC is a tuple (S,R,L) where: 

• S is a finite set of states; 

• R : S x S —> R>o is a transition rate matrix; 

• L : S — > 2 AP is a labeling function which associates each state with a set of atomic propositions. 

The transition rate matrix R assigns rates to each pair of states, which are used as parameters of the 
exponential distribution. A transition can occur between two states s and s' if R(s,s') > 0, and the 
probability of the transition being triggered within t time-units equals to 1 — e R ( s,s > 4 . If R(s,s') > for 
more than one state s', the first transition to be triggered determines the next state. Therefore, the choice 
of the successor state of s is probabilistic. The time spent in state s before any such transition occurs is 
exponentially distributed with E(s) = Y, s >esR{ s i s ')- Hence, the probability of moving to state s' is Ty , 
i.e., the probability that the delay of going from s to s' "finishes before" the delays of any other outgoing 
transition from s. A path in a CTMC is a sequence a in the form of ^o^o^i^i ' ■ ' with R(si,Sj + i) > and 
ti £ M>o for all i > 0. The amount of time spent in Sj is denoted by 

Corresponding to CTMC models, we use Continuous Stochastic Logic (CSL) to specify properties 
of built models. CSL, originally introduced by Aziz et al. O, provides a powerful means to specify both 
path-based and traditional state-based performance measures on CTMCs. 

Definition 2.2 The syntax of CSL is given as follows: 

::= true \ a | -.0 | A | P^p^U 1 ^} \ S^ p [^>] 
where a is an atomic proposition, ~G {<, <,>,>}, p € [0, 1], / is an interval o/'Mj.o cmd r,t E M>o- 
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CSL formulas are evaluated over the states of a CTMC. CSL includes the standard operators from propo- 
sitional logic: true (satisfied in all states); atomic propositions (a is true in states which are labelled with 
a); negation (-i<j) is true if <p is not); and conjunction (<j>i A 02 is true if both <p\ and 02 are true). Other 
standard boolean operators can be derived from these in the usual way. CSL also includes two proba- 
bilistic operators, P and S, both of which include a probability bound ~ p. A formula P^ p [y\ is true in a 
state s if the probability of the path formula y being satisfied from state s meets the bound ~ p. In this 
paper, we use a single type of path formula, F ! <p = trueU 1 ^, called an eventual formula, which is true 
for a path a if eventually becomes true for some time instant t £ I. Particularly, if the time interval is 
set to zero, e.g. F^Afy, the formula is true for a path a if becomes true at time instant t. The S operator 
is used to specify steady-state behavior of a CTMC. More precisely, S^ p [y} asserts that the steady-state 
probability of being in a state satisfying y meets the bound ~ p. 

2.2 The model checker PRISM 

PRISM ll26l is a model checking tool developed at the universities of Birmingham and Oxford. It al- 
lows one to model and analyze systems containing stochastic behaviors. PRISM supports three kinds of 
models: discrete-time Markov chains (DTMCs), continuous-time Markov chains (CTMCs) and Markov 
decision processes (MDPs). Analysis is performed through model checking such systems against prop- 
erties written in the probabilistic temporal logics PCTL if the model is a DTMC or an MDP, or CSL in 
the case of a CTMC, as well as their extensions for quantitative specifications and costs/rewards. 

In PRISM a model is composed of a number of modules that contain variables and can interact with 
each other. The values of the variables at any given time constitute the state of the module, and the local 
states of all modules decide the global state of the whole model. The behavior of a module, normally the 
changes in states which it can undergo, is specified by a set of guarded commands of the form: 

[a] g -)■ r : u; 

a is an action label in the style of process algebra, which introduces synchronization into the model. It 
can only be performed simultaneously by all modules that have an occurrence of action label a. If a 
transition does not have to synchronize with other transitions, then no action label needs to be provided 
for this transition. The symbol g is a predicate over all the variables in the system. A guarded command 
means that if the guard g is true, the system is updated according to u with rate r, which is corresponding 
to the transition rate of CTMC. A transition updates the value of variables by giving their new primed 
value with respect to their unprimed value. 

PRISM models can be augmented with information about rewards (or equivalently, costs). The tool 
can analyze properties which relate to the expected values of these rewards. A CTMC in PRISM can 
be augmented with two types of rewards: state reward associated with states which are accumulated 
in proportion to the time spent in the state, and transition reward associated with transitions which are 
accumulated each time the transition is taken. CSL is extended with quantitative costs/rewards as well, 
which is quite useful in analyzing the quantitative properties of a biological system, by introducing the 
R operator: 

R ::= R~ r [r l ] \ R~r[C*] \ R~ r [F<i>] I R~r[S] 

where ~G {<,<,>,>}, r,t G M>o and is a CSL formula. Intuitively, a state s satisfies R^ r [I =t ] if 
from s the expected state reward at time instance t meets the bound ~ r; a state s satisfies R^ r [C-'] if the 
expected reward accumulated up until t time units past satisfies ~ r; a state s satisfies R^ r [F(j>] if from s 
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the expected reward accumulated before a state satisfying (j) is reached meets the bound ~ r; and a state 
s satisfies /?^ r [5] if from s the long-run average expected reward satisfies ~ r. 

It is often useful to take a quantitative approach to probabilistic model checking, computing the 
actual probability that some behavior of a model is observed, rather than just verifying whether or not 
the probability is above or below a given bound. Hence, PRISM allows the P and S operators in CSL to 
take the following form: P=i[y] and S=?[i/a]. 

3 The PDGF Signaling Pathway 

Cell signaling is part of a complex system in cellular communication. It allows the cells to activate a 
large number of signaling molecules and to regulate their activity. In order to transfer a regulatory signal 
upon reception of a triggering stimulus, the signal is transformed into a chemical messenger within the 
signaling cell, e.g., via transfer of a phosphate group (phosphorylation) [22]. For further details on cell 
signaling see, for example, Il4ll22l. 

Platelet-Derived Growth Factor (PDGF), described approximately 30 years ago as a major mitogenic 
component of whole blood [37] is a growth factor that regulates cell growth and division. By binding 
to its receptor (PDGFR), it regulates many biological processes such as migration, survival and prolif- 
eration |fT9l . PDGFR is a receptor tyrosine kinase, which in general transfer upstream signals to many 
downstream signaling pathways by phosphorylation. Up to now, five PDGF ligands are known, PDGF- 
AA, -AB, -BB, -CC, -DD, interacting with three different types of PDGFR complexes, PDGFR-aa, -a/3 
and -j3/3. Each of the PDGFR subtypes has a different affinity to the different PDGF ligands IT331I . 

After PDGFRs couple with their respective ligands, phosphorylation of the receptor at specific tyro- 
sine residues will occur, thus enabling binding of signaling enzymes including Src, phosphatidylinositol 
3 kinase (PI3K), phospholipase Cy(PLCy) and SHP2 in the MAPK pathway at specific binding sites. The 
recruitment of these signaling enzymes to PDGFR is mediated via an intrinsic SH2 domain. The translo- 
cation of PI3K and PLC/ to the plasma membrane also increases their accessibility to their respective 
substrates. Moreover, recent findings suggest that PDGFR also has potential binding sites for CrkL II351L 
which will activate Rapl to positively influence c-Raf in the MAPK pathway [ 13], for Signal Transducer 
and Activator of Transcription (STAT), which might regulate the signal in parallel to the JAK-STAT path- 
way IT341 and also for cCbl, which promotes ubiquitination of PDGFR. cCbl is also considered to be one 
of the negative regulatory molecules in PDGF signal transduction [21 j. Based on intensive literature 
review, we have built a PDGF signal tranduction model in ODE format. This full model comprises 35 
molecules. The essential parts of the PDGF signaling pathway which are the coupling of PDGF ligand 
to PDGFR, the negative regulatory feedbacks on PDGFR and the activation of two main downstream sig- 
naling pathways, the MAPK and PI3K/Akt pathways, are extracted for analysis in PRISM (as shown in 
Figure [TJ. This simplified model contains only 17 molecules, and the analysis in this paper focuses on 
this simplified model. 

Figure [T] describes how signals are transduced in the PDGF pathway by activating or deactivating 
specific downstream pathways signaling molecules respectively. In this model, there are three inputs 
which are PDGFL (PDGF ligand), bPTEN, and bPDK. PDGFL is the node which represents the upstream 
molecule activating the whole network. PPX, and bPTEN are nodes which represent phosphatase enzymes 
in the cytoplasm that negatively regulate their targets. Lastly, bPDK, standing for basal activity of PDK, is 
the node that constantly gives a basal additional input to PDK node. This node is always active in order 
to activate the survival pathway to counteract apoptotic signal and keep the cell survive at a basal level. 
There are three different types of arrows in the network: blue arrows, green arrows and red arrows. The 
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Figure 1: The extracted PDGF signaling pathway (blue arrows: main pathway, green arrows: positive 
crosstalk, red arrows: negative regulatory) 



blue arrows represent the main activating interactions, indicating the two main downstream signaling 
pathways in the network. The MAPK pathway covers SHP2, Grb2S0S, GabSOS, Ras, c-Raf and MEK12. 
These molecules play a major role in the cellular proliferative circuit [ 16 ]. The PI3K/Akt pathway covers 
the molecules PI3K, PIP3, PDK, bPTEN, bPDK and Akt, also being important in viability circuit iTTol 
The green arrows represent positive crosstalk interactions to the molecules the arrows point to. Lastly, 
the red arrows represent either negative crosstalk interactions or other negative regulatory interactions. 
The molecules will become active after they have been activated by either blue or green arrows. In 
contrary, the molecules will become inactive after they have been deactivated by the red arrows or the 
basal phosphatase activities in the cell (not shown in the figure). 

PDGFR can be activated by PDGFL. The active PDGFR in turn activates three downstream molecules 
which are SHP2, PI3K and cCbl. Both SHP2 and cCbl assert a negative feedback to PDGF making it 
inactive. The three blue arrows connecting PDGFR to these downstream signaling molecules, so called 
mutant arrows, are the targets of system interventions both experimentally and computationally. The ex- 
perimental intervention can be performed by introducing a point mutation from tyrosine to phenylalanine 
at the specific recruitment site for the downstream signaling enzyme (Y720F for SHP2 recruitment site, 
YY731/742FF for PI3K recruitment site, and Y1018F for cCbl recruitment site), leading to the loss of 
signal capacity of the respective signaling pathway ll36l l3l 1291. Thus, the result of computational sim- 
ulation such as the relative activities at the steady state of downstream signaling molecules from these 
respective mutants can be validated experimentally in biological laboratories. 

Figure [2] contains the list of model reactions. Each node (molecule) is simplified to be in two states, 
either inactive or active (indicated by the suffix _act in the figure). All the reactions except for reactions 
9, 10 and 11 describe the molecules changing between the two states; while reactions 9, 10 and 11 
indicate the basal production, the basal degradation and the internalization following activation of PDGFR. 
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Figure 2: List of biochemical reactions in the PDGF signaling pathway ODE model 



For instance, reaction 2 describes that an active PPX gives a negative feedback to PDGFR, making it 
inactive. The parameters used in the model (all starting with letter k) in each reaction are the relative 
reaction rate constants normalized to the maximum value in each set of reactions. In this case, the 
values of parameters are ranging from to 1. The initial set of parameters' values used in this model 
are derived from literature (e.g., lTT2l [T] [371 HI HU) and our own experimental observations of PDGFRa 
mutants. Further parameter optimization based on experimental data is ongoing, especially focusing on 
the unknown crosstalk reaction strength. However, in this paper, the initial set of parameters has been 
used. This might not fully correlate to the actual biological system. Therefore, we focus in our analysis 
on the more qualitative aspects already conserved by the network structure. 



4 Modelling and Property Specifications in PRISM 
4.1 PRISM model 

We now describe how to build a PRISM model for the PDGF signaling pathway as presented in the 
previous section. Our model represents a single instance of the signaling pathway, meaning there can be 
at most one element of each molecule. In the single instance model, the molecules' steady state, which 
is expressed as a probability in PRISM, corresponds to the molecule density in the real-life experiments. 

Each of the nodes (or molecules) of the pathway except for the PDGFL in Figure [T] is represented 
by a separate PRISM module. Since PDGFL, bPTEN and bPDK remain the same after the reactions they 



are involved, we set them as const boolean values true. Figures 3(a) and 3(b) show the modules for 
PDGFR and PPX. In each of the modules, the status of the molecule is represented by a variable with 
the same name as the module. The variables can have values of either or 1 (PDGFR is an exception, 
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module PDGFR 

PDGFR : [0..2] init 0; 110 - inactive; 1 - active; 2 - degraded 

[] PDGFL & PDGFR=0 -> kon : (PDGFR '=1); 

[bkoffl] PDGFR=\ -> koffl : (PDGFR'=0); 

[bkoff2] PDGFR=l -> koffl : (PDGFR'=0); 

[bkubi] PDGFR=\ -> kubi : (PDGFR' =Q); 

[bkoffppx] PDGFR=\ -> koffppx : (PDGFR'=\); 

[bkl] PDGFR=\ -> kl : (PDGFR'=\); 



[MS] PDGFR=\ -> k5 : 
[bk6] PDGFR=\ -> k6 : 
[] PDGFR=0 -> kbasal 
[] PDGFR=2 -> kbasal 



( PDGFR =1); 
(PDGFR'=\); 

(PDGFR '=2); 1 1 PDGFR degraded 
(PDGFR '=0); 



[] PDGFR=\ -> kdeg : (PDGFR'=2); IIPDGFR_act degraded 



endmodule 



(a) PRISM module for PDGFR 



module 

PPAT: [0..1] init 1; 

[Mc#i] PPZ=7 (PPX'=l); 

[bkoffppx] PPX=1 (PPX'=0); 

endmodule 



rewards "PDGFRactive" 

PDGFR=\:\; 
endrewards 

(c) PRISM rewards 



(b) PRISM model for PPX 

Figure 3: PRISM modules and rewards 



since it can have value 2 since it can degrade), corresponding to the two states, inactive and active, of 
a molecule. Each command in the PRISM module represents a reaction in Figure [2] Interactions of 
multiple molecules are implemented by the synchronization between modules. More precisely, the same 
label is given to the commands which require synchronization in PRISM modules. For example, in 
Figures 3(a) and 3(b)| there are commands with label bkoffl in both of the modules PDGFR and PPX. 
The two commands are used to model the reaction (2) in Figure [2j which involves both PDGFR and PPX. 
It guarantees that the two commands (corresponding to one reaction) can only occur when both guards 
are satisfied. The reaction rate is assigned by the command in module PDGFR and hence the reaction rate 
of the command in module PPX is omitted. Totally we have modeled all the 17 molecules in 14 PRISM 
modules (PDGFL, bPTEN and bPDK are modeled as a constant). 



As mentioned in Section 2.2 PRISM models can be augmented with information about rewards. We 



construct rewards to calculate the time for a molecule being active. Figure 3(c) shows the rewards for 
calculating the active state of PDGFR. Each time PDGFR is in active state, one is added to the total time of 
PDGFR being active. Similarly, we build rewards structures for other molecules as well, including SHP2, 
Ras, MEK12, PIP3 and Akt. 



4.2 Property specifications 

There are three main goals for this study: (1) analyze the dynamics of PDGF induced signaling, (2) 
analyze the influence of the crosstalk reactions as defined in Section[3j and (3) analyze the importance of 
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Model 


States 


Transitions 


WildType 


589,824 


7,145,472 


SHP2Mutant 


36,864 


373,248 


PI3KMutant 


589,824 


7,096,320 


cCblMutant 


294,912 


3,357,696 


RasPI3KMutant 


589,824 


7,047,168 


Aktc-RafMutant 


786,423 


10,264,576 



Table 1: Model statistics in PRISM 



individual reaction on downstream signaling molecules. For the first goal, we study the signal transduc- 
tion properties of each mutant by removing the mutant arrows one by one and examine how the states of 
each molecule change accordingly at different time instances. We also examine the total time for each 
molecule being active. Moreover, it is interesting to study the activities of each molecule at the steady 
state as well. For the second goal, we do the comparison of probabilities for molecules to be active 
between different mutants by removing each of the crosstalk reactions. For the last one, we study how 
the steady state probabilties of molecule MEK12 and Akt change when a certain reaction is removed. 

Below we list properties of the PRISM model that we have analyzed to achieve our goals. Here, we 
use only the molecule PDGFR to illustrate the specification of the properties expressed as CSL formulas. 

• P =7 [FM PDGFR = 1] 

The probability that the molecule PDGFR is active at time instant t. 

. R^ DGFRacth ' e " } [C<=t] 

The expected time of PDGFR being active by time t . It refers to the reward structure "PDGFRactive" 



as defined in Figure 3(c) 



• 5 =? [ PDGFR = 1] 

The long-run probability that PDGFR is active. 



5 Results 



We use PRISM to construct the PDGF model described in Section 4.1 and analyze the set of properties 
listed in Section [4~2| As described in Section[3j 



For the first goal to analyze the dynamics of PDGF induced signaling, we first develop a base model 
representing the system, in which all the reactions in Figure [2] are included. Subsequently, we obtain 
the mutant models by removing the mutant arrows (as mentioned in Section [3]) one by one. More pre- 
cisely, we have developed four models in this experiment including the base model corresponding to the 
WildType condition. The second one, called SHP2Mutant, is obtained by removing the mutant arrow 
pointed to SHP2. The third and the fourth ones are PI3KMutant and cCblMutant. Following the same 
process, we remove the mutant arrows pointed to PI3K and cCbl separately in the last two models. The 
three removed mutant arrows correspond to the reactions (7), (8) and (6) in Figure [2| The size of the 
models in PRISM are shown in Table [T] (Table [T] also shows the size of models RasPI3KMutant and 
Aktc-RafMutant which will be discussed later). For each model, we compute the probability of each 
molecule being active at time instance t, which is summarized in Figure [4] 
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Figure|4]shows the probability of 6 out of the 17 molecules namely PDGFR, SHP2, Ras, PIP3, MEK12, 
and Akt. The 6 molecules are chosen according to their positions in the signaling pathway. We can see 



from Figure 4(a) that cCbl affects PDGFR more than the other two molecules (SHP2 and PI3K). This is 
due to strong negative feedback (red arrow) from cCbl to PDGFR in Figure [T] There is also a negative 
feedback from SHP2 to PDGFR. However, the reaction rate of the negative feedback from cCbl is 7 times 



as large as the one from SHP2, hence cCbl can affect PDGFR more than SHP2. The red curve in Figure 4(d) 
shows that not only PIP3 in PI3KMutant is less active than in other conditions, but also the time of it 
becoming active is delayed. In biological experiments, this delay is not observed. This is due to the 
fact that in the full PDGF signaling pathway model there is another molecule PLCg giving a fast positive 
crosstalk to PI3K, which is removed in our simplified PDGF model here. 

If we focus on the light blue curve of Figures 4(a)| - 4(f) we can see that cCbl has much impact on 
PDGFR but little impact on the other molecules. This is because after PDGFR activates SHP2 and PI3K, 
the states of the other molecules are determined by SHP2 and PI3K. Furthermore, from the green curve 
of Figures |4(b)| 4(c) and 4(e) we see that SHP2 has great effect on the status of the molecules on the 
MAPK pathway. However, the more downstream the molecule is, the less prominent the effect becomes 
(the effect to MEK12 is less than that to Ras). Due to the influence of the positive feedback from PI3K 
and PIP3, the molecules in the MAPK pathway can also become active without SHP2, which is the 



result from positive crosstalk interactions. Besides, the red curves of Figures 4(d) and 4(f) show that 
the influence of PI3K to the molecules in the PBK/Akt pathway, which is similar to that of SHP2 to the 
MAPK pathway. 

We analyze the long-run probability of the molecules being active after PDGFL stimulation (shown in 
Table [2]). The result demonstrates the activities of each molecule in this model at steady state. 



Molecule 


Probability 


Molecule 


Probability 


PDGFR 


0.22 


Grb2S0S 


0.55 


SHP2 


0.45 


Ras 


0.72 


GabSOS 


0.53 


MEK12 


0.77 


c-Raf 


0.63 


PIP3 


0.72 


PI3K 


0.62 


Akt 


0.82 


PDK 


0.83 


MT0R 


0.84 


cCbl 


0.47 


PPX 


0.00 



Table 2: Steady state probabilities of molecules 



PRISM supports reward properties (see Section|4]). Figure[5]shows the expected time of six molecules 
being active by time instant t. These six molecules are the same as in Figure [4j All the six curves tend to 
be linear after time instant 12, which shows that the state of the 6 molecules start to be in a steady state 
after time instant 12. 

For the second goal to analyze the influence of the crosstalk reactions, the experiment is also per- 
formed as in silico genetics. After developing the base model, we get the model variants by removing 
one crosstalk arrow at a time. Totally, there are four positive crosstalk reactions (green arrows) and two 
negative crosstalk reactions (red arrows pointing from Akt). In this analysis, we focus on the positive 
crosstalk from Ras to PI3K and the negative crosstalk from Akt to c-Raf. More precisely, we develop 
three models in this experiment. The first one is the base WildType model with all reactions included. 
The second one is the RasPI3KMutant model, in which the positive feedback from Ras to PI3K is re- 
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Figure 5: Expected time of being active by time t 




moved. The last one is the Aktc-Raf Mutant model, in which the negative feedback from Akt to c-Raf 
is removed. After building the three models, we analyze the influences of the two crosstalks by com- 
paring the WildType model to the two mutant models, respectively. The size of each PRISM model is 
summarized in Table [T] 

Figure [6] shows the results of the comparison. In both sub-figures, we compare the molecules which 
are directly related to the crosstalk arrows and the molecules at the end of the related main signaling 
pathways. In Figure 6(a) we compare in each model the probabilities for molecules PI3K and Akt 
being active, which is at the end of the PI3K/Akt pathway; while in Figure 6(b) we have chosen the 
molecules c-Raf and MEK12, which is at the end of the MAPK pathway. The dark blue and green curves 
in Figure [6(a) show how the green arrow influences the status of PI3K. If there is no positive crosstalk, 
the probability of PI3K being active (green curve) becomes smaller. The two curves are coincident for 
a time period because the molecule Ras needs some time to become active before it can activate PI3K. 
The red and light blue curves show that the influence of the positive feedback to Akt is smaller than that 
to PI3K. Just like the situation in Figures 14(b) 1 14(c)! and |4(e)l the more downstream the molecules are, 



Yuan et al. 



77 



the smaller the influence becomes. The results in Figure 6(b) are similar: the dark blue and green curves 
show that the influence the negative crosstalk gives to Ras is significant; while the positive crosstalk has 
little influence on MEK12. 

For the last goal to analyze the importance of individual reaction on downstream signaling molecules, 
we compute the steady state probabilities of Akt and MEK12 in 31 different models, each of which is 
obtained by removing one reaction from the wildtype model. For example, 'PIP3-GabSOS' model is 
obtained by removing reaction No. 20 in Figure [2]). The results are shown in Figure [7] We draw two 
lines in the figure to divide the axis into four areas, letting the dot of the wildtype model lie at the 
cross of the two lines. Apparently, the dots in different areas show that the removed reactions can bring 
different influence to the steady state probabilities of Akt and MEK12. The reactions (corresponding to 
the dots) in area 1 can decrease the steady state probabilities of both Akt and MEK12, while those in area 
3 can increase both probabilities. The reactions in area 2 can decrease the steady state probability of 
Akt and increase the one of MEK12; while reactions grouped in area 4 lead to the opposite effect. For 
those dots lying on the horizontal line, the corresponding removed reactions have little impact on the 
steady state probability of Akt. Similarly, the reactions on the vertical line have little impact on MEK12. 
These observations have biological implications. Akt and MEK12 are downstream molecules in the signal 
transduction process that regulate different cellular functions. Signal from Akt keeps the cells to survive 
from apoptosis and signal form MEK12 regulates the cells growth and proliferation. In cancer, both of 
these two main pathways (see Figure [T} are more active so they drive the cells to keep growing and 
dividing in an uncontrolled manner. Therefore, if we could find the targets to control these two signaling 
molecules to be at a desirable level, it would be beneficial for cancer therapeutic development. 



6 Conclusion and Future Work 

In this paper, we have given a detailed description of using the probabilistic model checker PRISM to 
study the PDGF signaling pathway. Based on intensive literature reviews, we have built a PDGF signal 
transduction model in ODE format and extracted it for the analysis in a probabilistic model checker. 
We focused on a simplified model of the PDGF signaling pathway which contains 17 molecules while 
the full model has 35 molecules. We have analyzed the dynamics of PDGF induced signaling, the 
influence of crosstalk reactions in the signaling pathway and the importance of individual reaction on 
downstream signaling molecules. Our experiment results show that quantitative verification can provide 
us a better understanding of the PDGF signaling pathway, especially the result discussed in the end of 
Sect. [5] potentially can give rise to better behavior prediction of the pathway. (The PRISM model and 
property specifications can be found at satoss . uni . lu/jun/models/PDGF . zip.) 

Except for model checking, we have also analyzed the model using ODE simulations. The results for 
the two molecules, PDGFR and Ras, are shown in Figure 8(a)| and Figure 8(b) respectively. We can see 



that probabilistic verification and ODE simulations produce comparable results: the shape of the curves 
are similar while the maximum probabilities and the steady state probabilities differ a little. The reason 
of the difference might be that the approach based on verification leads to more accurate results when 
the number of molecules is small as indicated in [27]. In general, ODE analysis perform well for larger 
number of molecules, but not for small numbers, as time trajectories for average concentrations can be 
misleading if numbers of molecules small. We plan to have a detailed comparison between our analysis 
using probabilisitic verification techniques and the analysis using ODE simulations. As indicated in 
the end of Section [3] our model is not intended to be fully accurate. The reaction rates, especially 
the crosstalk reaction rates, are based on literature reviewing and are only relative values. Currently, 
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Figure 7: Steady state probabilities of Akt and MEK12 in different models 
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experiments in a biological laboratory are performed to get more precise values for these reaction rates. 
We hope techniques like parametric verification for Markov chains (e.g., |[T5l[T4l ) can help to synthesize 
values which are consistent with the lab experimental results. 
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